#### Does mislabeling COVID-19 elicit the perception of threat and reduce blame?
#### Chengxin Xu & Yixin Liu
#### JBPA Replication file: appendix tables and figures

# install.packages(c('list','ggplot2','plotly','ggpubr','egg','simpleboot','cowplot','Rmisc','list','xtable','stargazer'))

library(list)
library(ggplot2)
library(plotly)
library(ggpubr)
library(egg)
library(simpleboot)
library(cowplot)
library(Rmisc)
library(list)
library(xtable)
library(stargazer)

##import data
setwd()
exp<-read.csv("JBPA_Xu&Liu.csv")

####Subgroup Datasets####
exp.con<-subset(exp,con==1) #subgroup: conservatives
exp.lib<-subset(exp,lib==1) #subgroup: liberals
exp.mod<-subset(exp,mod==1) #subgroup: moderates
exp.dem<-subset(exp,Dem==1) #subgroup: Democrats
exp.rep<-subset(exp,Rep==1) #subgroup: Republicans
exp.ind<-subset(exp,Ind==1) #subgroup: Independents
exp.hs<-subset(exp,hsocial==1) #subgroup: High media use
exp.ls<-subset(exp,hsocial==0) #subgroup: low media use
exp.covid.c<-subset(exp,covid==0) #subgroup: covid group 3 item
exp.covid.t<-subset(exp,covid==1) #subgroup: covid group 4 item
exp.virus.c<-subset(exp,virus==0) #subgroup: virus group 3 item
exp.virus.t<-subset(exp,virus==1) #subgroup: virus group 4 item

####Figure Aesthetic####
myShape <- c(19,15,17)
pos2 <- position_dodge(width=0.2)
pos3 <- position_dodge(width=0.5)
limits <- aes(ymax = upr, ymin = lwr)

####Appendix A####
## Characteristics of sample

##randomization check
sex.aov <- aov(sex ~ group, data = exp)
summary(sex.aov)
male<-summary(sex.aov)[[1]][,5][1]
female.aov <- aov(female ~ group, data = exp)
summary(female.aov)
female<-summary(female.aov)[[1]][,5][1]
young.aov <- aov(young ~ group, data = exp)
summary(young.aov)
young<-summary(young.aov)[[1]][,5][1]
mid.aov <- aov(midage ~ group, data = exp)
summary(mid.aov)
mid<-summary(mid.aov)[[1]][,5][1]
old.aov <- aov(old ~ group, data = exp)
summary(old.aov)
old<-summary(old.aov)[[1]][,5][1]
wht.aov <- aov(white ~ group, data = exp)
summary(wht.aov)
wht<-summary(wht.aov)[[1]][,5][1]
blk.aov <- aov(blk ~ group, data = exp)
summary(blk.aov)
blk<-summary(blk.aov)[[1]][,5][1]
his.aov <- aov(his ~ group, data = exp)
summary(his.aov)
his<-summary(his.aov)[[1]][,5][1]
asian.aov <- aov(asian ~ group, data = exp)
summary(asian.aov)
asi<-summary(asian.aov)[[1]][,5][1]
other.aov <- aov(other ~ group, data = exp)
summary(other.aov)
oth<-summary(other.aov)[[1]][,5][1]
hwor.aov <- aov(hwor ~ group, data = exp)
summary(hwor.aov)
hwor<-summary(hwor.aov)[[1]][,5][1]
lwor.aov <- aov(lwor ~ group, data = exp)
summary(lwor.aov)
lwor<-summary(lwor.aov)[[1]][,5][1]
hsocial.aov <- aov(hsocial ~ group, data = exp)
summary(hsocial.aov) ##*
hso<-summary(hsocial.aov)[[1]][,5][1]
lsocial.aov <- aov(lsocial ~ group, data = exp)
summary(lsocial.aov) ##*
lso<-summary(lsocial.aov)[[1]][,5][1]
dem.aov <- aov(Dem ~ group, data = exp)
summary(dem.aov)
dem<-summary(dem.aov)[[1]][,5][1]
rep.aov <- aov(Rep ~ group, data = exp)
summary(rep.aov)
rep<-summary(rep.aov)[[1]][,5][1]
ind.aov <- aov(Ind ~ group, data = exp)
summary(ind.aov) 
ind<-summary(ind.aov)[[1]][,5][1]
con.aov <- aov(con ~ group, data = exp)
summary(con.aov)
con<-summary(con.aov)[[1]][,5][1]
lib.aov <- aov(lib ~ group, data = exp)
summary(lib.aov) 
lib<-summary(lib.aov)[[1]][,5][1]
mod.aov <- aov(mod ~ group, data = exp)
summary(mod.aov)
mod<-summary(mod.aov)[[1]][,5][1]
low.aov <- aov(low ~ group, data = exp)
summary(low.aov)
low<-summary(low.aov)[[1]][,5][1]
med.aov <- aov(med ~ group, data = exp)
summary(med.aov)
med<-summary(med.aov)[[1]][,5][1]
high.aov <- aov(high ~ group, data = exp)
summary(high.aov)
high<-summary(high.aov)[[1]][,5][1]

####Table A.1####
des_stat<-matrix(NA,23,11)
rownames(des_stat)<-c("Male","Female","Age: 18-29","Age: 30-49",
                      "Age: 50 and older","White","Black","Hispanic",
                      "Asian","Other","COVID-19: worried","COVID-19: not worried",
                      "More than half time Twitter/Facebook",
                      "Less than half time Twitter/Facebook",
                      "Democrat","Republican","Independent",
                      "liberal","Conservative","Moderate",
                      "Income: Less than $25,000",
                      "Income: $25,000 to $74,999",
                      "Income: $75,000 or more")
colnames(des_stat)<-c("Frequency","%",
                      "Frequency","%",
                      "Frequency","%",
                      "Frequency","%",
                      "Frequency","%","Pr(>F)")
tol.sum<-round(apply(exp[,c(41:63)],2,sum),digits=0)
cc.sum<-round(apply(exp.covid.c[,c(41:63)],2,sum),digits=0)
ct.sum<-round(apply(exp.covid.t[,c(41:63)],2,sum),digits=0)
tc.sum<-round(apply(exp.virus.c[,c(41:63)],2,sum),digits=0)
tt.sum<-round(apply(exp.virus.t[,c(41:63)],2,sum),digits=0)

# percentage function
pcentFun <- function(x) {
  res <- x
  100 * (sum(res) / length(res))
}

tol.p<-round(apply(exp[,c(41:63)],2,pcentFun),2)
cc.p<-round(apply(exp.covid.c[,c(41:63)],2,pcentFun),2)
ct.p<-round(apply(exp.covid.t[,c(41:63)],2,pcentFun),2)
tc.p<-round(apply(exp.virus.c[,c(41:63)],2,pcentFun),2)
tt.p<-round(apply(exp.virus.t[,c(41:63)],2,pcentFun),2)
des_stat[,1]<-tol.sum
des_stat[,3]<-cc.sum
des_stat[,5]<-ct.sum
des_stat[,7]<-tc.sum
des_stat[,9]<-tt.sum
des_stat[,2]<-tol.p
des_stat[,4]<-cc.p
des_stat[,6]<-ct.p
des_stat[,8]<-tc.p
des_stat[,10]<-tt.p
des_stat[,11]<-c(male,female,young,mid,old,wht,blk,his,asi,oth,hwor,lwor,
                 hso,lso,dem,rep,ind,lib,con,mod,low,med,high)
des_stat[,11]<-round(des_stat[,11],digits = 2)
print(xtable(x=des_stat))
write.table(des_stat, file = "a1.txt", sep = ",", quote = FALSE, row.names = T)

####Appendix B####
##Design effect tests for list experiment

test.value <- ict.test(exp$ict, exp$list, 
                       alpha = 0.05, 
                       J = 3, gms = TRUE)
print(test.value)

test.value.r <- ict.test(exp.rep$ict, exp.rep$list, 
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.r)

test.value.d <- ict.test(exp.dem$ict, exp.dem$list, 
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.d)

test.value.i <- ict.test(exp.ind$ict, exp.ind$list, 
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.i)

test.value.c <- ict.test(exp.con$ict, exp.con$list,
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.c)

test.value.l <- ict.test(exp.lib$ict, exp.lib$list,
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.l)

test.value.m <- ict.test(exp.mod$ict, exp.mod$list,
                         alpha = 0.05, 
                         J = 3, gms = TRUE)
print(test.value.m)

####Table B.1####
design.test<-rbind(test.value$pi.table,c(test.value$p,NA))
rownames(design.test)[rownames(design.test)==""]<-"Bonferroni-corrected p-value"
design.test<-round(design.test,2)
design.test
print(xtable(x=design.test))
write.table(design.test, file = "b1.txt", sep = ",", quote = FALSE, row.names = T)

####Table B.2####
design.test.party<-cbind(test.value.d$pi.table,
                         test.value.i$pi.table,
                         test.value.r$pi.table)
design.test.party.p<-cbind(test.value.d$p,NA,
                           test.value.i$p,NA,
                           test.value.r$p,NA)
design.test.party<-rbind(design.test.party,design.test.party.p)
rownames(design.test.party)[rownames(design.test.party)=="Sensitive Item 1"]<-"Bonferroni-corrected p-value"
design.test.party<-round(design.test.party,2)
design.test.party
print(xtable(x=design.test.party))
write.table(design.test.party, file = "b2.txt", sep = ",", quote = FALSE, row.names = T)

####Table B.3####
design.test.ideo<-cbind(test.value.c$pi.table,
                        test.value.m$pi.table,
                        test.value.l$pi.table)
design.test.ideo.p<-cbind(test.value.c$p,NA,
                          test.value.m$p,NA,
                          test.value.l$p,NA)
design.test.ideo<-rbind(design.test.ideo,design.test.ideo.p)
rownames(design.test.ideo)[rownames(design.test.ideo)=="Sensitive Item 1"]<-"Bonferroni-corrected p-value"
design.test.ideo<-round(design.test.ideo,2)
design.test.ideo
print(xtable(x=design.test.ideo))
write.table(design.test.ideo, file = "b3.txt", sep = ",", quote = FALSE, row.names = T)

####Appendix C####
##The Heterogeneity of the “Chinese Virus” Label Treatment Effect 
###on Social Media Use

# high social media use
nls.hs <- ictreg(ict~ch,
                 treat = "list",
                 data=exp.hs,
                 J=3, method = "nls")

summary(nls.hs)

exp.virus.hs <- exp.covid.hs <- exp.hs
exp.virus.hs[, "ch"] <- 1
exp.covid.hs[, "ch"] <- 0

avg.pred.diff.hs <- predict(nls.hs,
                            newdata = exp.virus.hs, 
                            newdata.diff = exp.covid.hs,
                            se.fit = TRUE, avg = TRUE, 
                            interval="confidence",
                            level=0.90)

avg.pred.diff.hs$fit
covid.hs.est<-avg.pred.diff.hs$fit[2,1]
virus.hs.est<-avg.pred.diff.hs$fit[1,1]
diff.hs.est<-avg.pred.diff.hs$fit[3,1]
covid.hs.upr<-avg.pred.diff.hs$fit[2,3]
covid.hs.lwr<-avg.pred.diff.hs$fit[2,2]
virus.hs.upr<-avg.pred.diff.hs$fit[1,3]
virus.hs.lwr<-avg.pred.diff.hs$fit[1,2]
diff.hs.upr<-avg.pred.diff.hs$fit[3,3]
diff.hs.lwr<-avg.pred.diff.hs$fit[3,2]

# low social media use
nls.ls <- ictreg(ict~ch,
                 treat = "list",
                 data=exp.ls,
                 J=3, method = "nls")

summary(nls.ls)

exp.virus.ls <- exp.covid.ls <- exp.ls
exp.virus.ls[, "ch"] <- 1
exp.covid.ls[, "ch"] <- 0

avg.pred.diff.ls <- predict(nls.ls,
                            newdata = exp.virus.ls, 
                            newdata.diff = exp.covid.ls,
                            se.fit = TRUE, avg = TRUE, 
                            interval="confidence",
                            level=0.90)

avg.pred.diff.ls$fit
covid.ls.est<-avg.pred.diff.ls$fit[2,1]
virus.ls.est<-avg.pred.diff.ls$fit[1,1]
diff.ls.est<-avg.pred.diff.ls$fit[3,1]
covid.ls.upr<-avg.pred.diff.ls$fit[2,3]
covid.ls.lwr<-avg.pred.diff.ls$fit[2,2]
virus.ls.upr<-avg.pred.diff.ls$fit[1,3]
virus.ls.lwr<-avg.pred.diff.ls$fit[1,2]
diff.ls.upr<-avg.pred.diff.ls$fit[3,3]
diff.ls.lwr<-avg.pred.diff.ls$fit[3,2]

####Figure C.1####
media.df<-data.frame(est=c(covid.hs.est,virus.hs.est,diff.hs.est,
                           covid.ls.est,virus.ls.est,diff.ls.est),
                     upr=c(covid.hs.upr,virus.hs.upr,diff.hs.upr,
                           covid.ls.upr,virus.ls.upr,diff.ls.upr),
                     lwr=c(covid.hs.lwr,virus.hs.lwr,diff.hs.lwr,
                           covid.ls.lwr,virus.ls.lwr,diff.ls.lwr),
                     mtd=c("a","b","c",
                           "a","b","c"),
                     trt=c("High social media useage","High social media useage",
                           "High social media useage",
                           "Low social media useage","Low social media useage",
                           "Low social media useage"))
media.df$trt <- factor(media.df$trt, 
                       levels=c("High social media useage","Low social media useage"))

fc1.plot <- ggplot(media.df, 
                   aes(y=est, x=trt, group=mtd, shape=mtd))+
  theme_article()+geom_hline(yintercept = 0)+
  scale_shape_manual(labels=c("COVID-19", "Chinese Virus", "Difference"),
                     values =myShape)
fc1.plot <- fc1.plot + geom_point(size=4, position=pos3) + 
  geom_errorbar(limits, width=0, colour="black", position = pos3) + 
  labs(x=element_blank(),y="Estimated Proportion/Difference") + 
  theme(legend.title=element_blank(),axis.text=element_text(size=10),
        legend.justification=c(0,0), legend.position=c(0,0))

fc1.plot #Figure C.1
# Export Figure C.1
ggsave(file = "fc1.eps", width = 6, height = 3, dpi = 666)

####Appendix D####
##Overt Perception of Threat of Chinese Immigrants

direct<-glm(direct~ch,data=exp,family = "binomial")
summary(direct)
direct.l<-glm(direct~ch,data=exp.lib,family = "binomial")
summary(direct.l)
direct.c<-glm(direct~ch,data=exp.con,family = "binomial")
summary(direct.c)
direct.m<-glm(direct~ch,data=exp.mod,family = "binomial")
summary(direct.m)

####Table D.1####
stargazer(direct,direct.l,direct.c,direct.m,
          style = "apsr",
          covariate.labels=c("Chinese Virus"),
          column.labels=c("Overall",
                          "Liberal","Conservative","Moderate"),
          dep.var.labels.include = FALSE,
          intercept.bottom=TRUE,
          align=TRUE,no.space=TRUE)

####Appendix E####
##Robustness check drop manipulation check failure sample

exp<-subset(exp,exp$mcr==exp$ch)
exp.con<-subset(exp,con==1) #subgroup: conservatives
exp.lib<-subset(exp,lib==1) #subgroup: liberals
exp.mod<-subset(exp,mod==1) #subgroup: moderates

# Figure E.1 upper panel (item count)
ict.df<-summarySE(exp, measurevar="ict", 
                  groupvars=c("Item","ch1"),
                  conf.interval = 0.90)

# plot E.1.1
fe1.1.plot <- ggplot(ict.df, 
                    aes(y=ict, x=ch1, group=Item, shape=Item))+
  scale_shape_manual(values =c(1,16))+
  scale_x_discrete(labels=c("COVID-19","Chinese Virus"))+
  theme_article()
fe1.1.plot <- fe1.1.plot + geom_point(size=4, position=pos2) + 
  geom_errorbar(aes(ymin=ict-ci, ymax=ict+ci), width=0, colour="black", 
                position = pos2) + 
  labs(x=element_blank(),y="List Item Count") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(1,1), legend.position=c(1,1))+
  guides(colour=FALSE, shape = guide_legend(reverse=T))
fe1.1.plot

# Figure E.1 lower panel (proportion estimation)
# List estimation
st.mod.c = t.test(exp$ict[exp$group=="covid4"],
                  exp$ict[exp$group=="covid3"],
                  conf.level=0.90)
stc.est<-st.mod.c$estimate[1]-st.mod.c$estimate[2]
upper.stc<-st.mod.c$conf.int[2]
lower.stc<-st.mod.c$conf.int[1]
stc.ci<-upper.stc-stc.est

st.mod.t = t.test(exp$ict[exp$group=="virus4"],
                  exp$ict[exp$group=="virus3"],
                  conf.level=0.90)
stt.est<-st.mod.t$estimate[1]-st.mod.t$estimate[2]
upper.stt<-st.mod.t$conf.int[2]
lower.stt<-st.mod.t$conf.int[1]
stt.ci<-upper.stt-stt.est

# Direct estimation
exp.direct<-subset(exp,exp$direct==1|exp$direct==0)
direct.df<-summarySE(exp.direct, measurevar="direct", 
                     groupvars=c("ch1"),conf.interval = 0.90)

st.df<-data.frame(ch1=c("a","b"),N=c(NA,NA),direct=c(stc.est,stt.est),
                  sd=c(NA,NA),se=c(NA,NA),ci=c(stc.ci,stt.ci))
ict2.df<-rbind(direct.df,st.df)
ict2.df$mode<-c("Direct","Direct","List","List")

# plot E.1.2
fe1.2.plot <- ggplot(ict2.df, 
                    aes(y=direct, x=ch1, group=mode, shape=mode))+
  scale_shape_manual(values =c(2,17))+
  scale_x_discrete(labels=c("COVID-19","Chinese Virus"))+
  theme_article()+geom_hline(yintercept = 0)
fe1.2.plot <- fe1.2.plot + geom_point(size=4, position=pos2) + 
  geom_errorbar(aes(ymin=direct-ci, ymax=direct+ci), width=0, colour="black", 
                position = pos2) + 
  labs(x=element_blank(),y="Estimated Proportion") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(1,1), legend.position=c(1,1))+
  guides(colour = FALSE, shape = guide_legend(reverse=T))
fe1.2.plot

####Figure E.1####
##List Item Count and Proportion 
##of Participants Who Reported Perceived Threat
fe1.plot<-ggarrange(fe1.1.plot+theme(axis.text.x = element_blank(),
                                   axis.ticks.x = element_blank()),
                   fe1.2.plot,
                   ncol = 1,nrow = 2)

fe1.plot<-plot_grid(fe1.plot)
fe1.plot #Figure E.1
# export Figure E.1
ggsave(file = "fe1.eps", width = 4.5, height = 3, dpi = 666)

####Figure E.2####
##The “Chinese Virus” Label Treatment Effect

# DIM (bootstrap CI)
exp.covid<-subset(exp,ch==0)
exp.virus<-subset(exp,ch==1)

set.seed(123)
bootstrapped.c <- two.boot(exp.covid$ict.c[exp.covid$covid==1], 
                           exp.covid$ict.c[exp.covid$covid==0],
                           mean,
                           10000)
boot.mean_diff.covid <- data.frame(bootstrapped.c$t)
colnames(boot.mean_diff.covid) <- 'mean_diffs_c'

bootstrapped.t <- two.boot(exp.virus$ict.t[exp.virus$virus==1], 
                           exp.virus$ict.t[exp.virus$virus==0],
                           mean,
                           10000)

boot.mean_diff.virus <- data.frame(bootstrapped.t$t)
colnames(boot.mean_diff.virus) <- 'mean_diffs_v'
diffindiff <- boot.mean_diff.virus$mean_diffs_v-boot.mean_diff.covid$mean_diffs_c
mean.diffindiff <- mean(diffindiff)
lower.diffindiff <- quantile(diffindiff, probs = 0.05)
upper.diffindiff <- quantile(diffindiff, probs = 0.95)
cbind(lower.diffindiff, mean.diffindiff, upper.diffindiff)

# MLE
ml.results <- ictreg(ict~ch,
                     treat = "list",
                     data=exp,
                     J=3, method = "ml")

summary(ml.results)

exp.virus <- exp.covid <- exp
exp.virus[, "ch"] <- 1
exp.covid[, "ch"] <- 0

avg.pred.diff.mle <- predict(ml.results,
                             newdata = exp.virus, 
                             newdata.diff = exp.covid,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)
avg.pred.diff.mle$fit
covid.mle.est<-avg.pred.diff.mle$fit[2,1]
virus.mle.est<-avg.pred.diff.mle$fit[1,1]
diff.mle.est<-avg.pred.diff.mle$fit[3,1]
covid.mle.upr<-avg.pred.diff.mle$fit[2,3]
covid.mle.lwr<-avg.pred.diff.mle$fit[2,2]
virus.mle.upr<-avg.pred.diff.mle$fit[1,3]
virus.mle.lwr<-avg.pred.diff.mle$fit[1,2]
diff.mle.upr<-avg.pred.diff.mle$fit[3,3]
diff.mle.lwr<-avg.pred.diff.mle$fit[3,2]

# NLS
nls.results <- ictreg(ict~ch,
                      treat = "list",
                      data=exp,
                      J=3, method = "nls")

summary(nls.results)

exp.virus <- exp.covid <- exp
exp.virus[, "ch"] <- 1
exp.covid[, "ch"] <- 0

avg.pred.diff.nls <- predict(nls.results,
                             newdata = exp.virus, 
                             newdata.diff = exp.covid,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

avg.pred.diff.nls$fit
covid.nls.est<-avg.pred.diff.nls$fit[2,1]
virus.nls.est<-avg.pred.diff.nls$fit[1,1]
diff.nls.est<-avg.pred.diff.nls$fit[3,1]
covid.nls.upr<-avg.pred.diff.nls$fit[2,3]
covid.nls.lwr<-avg.pred.diff.nls$fit[2,2]
virus.nls.upr<-avg.pred.diff.nls$fit[1,3]
virus.nls.lwr<-avg.pred.diff.nls$fit[1,2]
diff.nls.upr<-avg.pred.diff.nls$fit[3,3]
diff.nls.lwr<-avg.pred.diff.nls$fit[3,2]

## Combine bootmean, MLE, NLS plots [MAIN]
boot.mle.nls<-data.frame(est=c(stc.est,covid.mle.est,covid.nls.est,
                               stt.est,virus.mle.est,virus.nls.est,
                               mean.diffindiff,diff.mle.est,diff.nls.est),
                         upr=c(upper.stc,covid.mle.upr,covid.nls.upr,
                               upper.stt,virus.mle.upr,virus.nls.upr,
                               upper.diffindiff,diff.mle.upr,diff.nls.upr),
                         lwr=c(lower.stc,covid.mle.lwr,covid.nls.lwr,
                               lower.stt,virus.mle.lwr,virus.nls.lwr,
                               lower.diffindiff,diff.mle.lwr,diff.nls.lwr),
                         mtd=c("DIM","MLE","NLS","DIM","MLE","NLS","DIM","MLE","NLS"),
                         trt=c("COVID-19","COVID-19","COVID-19",
                               "Chinese Virus","Chinese Virus","Chinese Virus",
                               "Difference","Difference","Difference"))
boot.mle.nls$trt <- factor(boot.mle.nls$trt, levels=c("COVID-19","Chinese Virus","Difference"))

fe2.plot <- ggplot(boot.mle.nls, 
                  aes(y=est, x=trt, group=mtd, shape=mtd))+
  scale_shape_manual(labels=c("DIM", "MLE", "NLS"),
                     values =myShape)+
  theme_article()+geom_hline(yintercept = 0)
fe2.plot <- fe2.plot + geom_point(size=4, position=pos3) + 
  geom_errorbar(limits, width=0, colour="black", position = pos3) + 
  labs(x=element_blank(),y="Estimated Proportion/Difference") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(0,0), legend.position=c(0,0))

fe2.plot # Figure E.2
# Export Figure E.2
ggsave(file = "fe2.eps", width = 6, height = 3, dpi = 666)

####Figure E.3####
# Liberal
nls.lib <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.lib,
                  J=3, method = "nls")

summary(nls.lib)

exp.virus.l <- exp.covid.l <- exp.lib
exp.virus.l[, "ch"] <- 1
exp.covid.l[, "ch"] <- 0

avg.pred.diff.lib <- predict(nls.lib,
                             newdata = exp.virus.l, 
                             newdata.diff = exp.covid.l,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

avg.pred.diff.lib$fit
avg.pred.diff.lib$se.fit
covid.lib.est<-avg.pred.diff.lib$fit[2,1]
virus.lib.est<-avg.pred.diff.lib$fit[1,1]
diff.lib.est<-avg.pred.diff.lib$fit[3,1]
covid.lib.upr<-avg.pred.diff.lib$fit[2,3]
covid.lib.lwr<-avg.pred.diff.lib$fit[2,2]
virus.lib.upr<-avg.pred.diff.lib$fit[1,3]
virus.lib.lwr<-avg.pred.diff.lib$fit[1,2]
diff.lib.upr<-avg.pred.diff.lib$fit[3,3]
diff.lib.lwr<-avg.pred.diff.lib$fit[3,2]

# Moderate
nls.mod <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.mod,
                  J=3, method = "nls")

summary(nls.mod)

exp.virus.m <- exp.covid.m <- exp.mod
exp.virus.m[, "ch"] <- 1
exp.covid.m[, "ch"] <- 0

avg.pred.diff.mod <- predict(nls.mod,
                             newdata = exp.virus.m, 
                             newdata.diff = exp.covid.m,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

avg.pred.diff.mod$fit
covid.mod.est<-avg.pred.diff.mod$fit[2,1]
virus.mod.est<-avg.pred.diff.mod$fit[1,1]
diff.mod.est<-avg.pred.diff.mod$fit[3,1]
covid.mod.upr<-avg.pred.diff.mod$fit[2,3]
covid.mod.lwr<-avg.pred.diff.mod$fit[2,2]
virus.mod.upr<-avg.pred.diff.mod$fit[1,3]
virus.mod.lwr<-avg.pred.diff.mod$fit[1,2]
diff.mod.upr<-avg.pred.diff.mod$fit[3,3]
diff.mod.lwr<-avg.pred.diff.mod$fit[3,2]

# Conservative
nls.con <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.con,
                  J=3, method = "nls")

summary(nls.con)

exp.virus.c <- exp.covid.c <- exp.con
exp.virus.c[, "ch"] <- 1
exp.covid.c[, "ch"] <- 0

avg.pred.diff.con <- predict(nls.con,
                             newdata = exp.virus.c, 
                             newdata.diff = exp.covid.c,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

avg.pred.diff.con$fit
covid.con.est<-avg.pred.diff.con$fit[2,1]
virus.con.est<-avg.pred.diff.con$fit[1,1]
diff.con.est<-avg.pred.diff.con$fit[3,1]
covid.con.upr<-avg.pred.diff.con$fit[2,3]
covid.con.lwr<-avg.pred.diff.con$fit[2,2]
virus.con.upr<-avg.pred.diff.con$fit[1,3]
virus.con.lwr<-avg.pred.diff.con$fit[1,2]
diff.con.upr<-avg.pred.diff.con$fit[3,3]
diff.con.lwr<-avg.pred.diff.con$fit[3,2]

## Combine ideology plots
ideo.df<-data.frame(est=c(covid.lib.est,virus.lib.est,diff.lib.est,
                          covid.mod.est,virus.mod.est,diff.mod.est,
                          covid.con.est,virus.con.est,diff.con.est),
                    upr=c(covid.lib.upr,virus.lib.upr,diff.lib.upr,
                          covid.mod.upr,virus.mod.upr,diff.mod.upr,
                          covid.con.upr,virus.con.upr,diff.con.upr),
                    lwr=c(covid.lib.lwr,virus.lib.lwr,diff.lib.lwr,
                          covid.mod.lwr,virus.mod.lwr,diff.mod.lwr,
                          covid.con.lwr,virus.con.lwr,diff.con.lwr),
                    mtd=c("a","b","c",
                          "a","b","c",
                          "a","b","c"),
                    trt=c("Liberal","Liberal","Liberal",
                          "Moderate","Moderate","Moderate",
                          "Conservative","Conservative","Conservative"))
ideo.df$trt <- factor(ideo.df$trt, 
                      levels=c("Liberal","Moderate","Conservative"))

fe3.plot <- ggplot(ideo.df, 
                  aes(y=est, x=trt, group=mtd, shape=mtd))+
  scale_shape_manual(labels=c("COVID-19", "Chinese Virus", "Difference"),
                     values =myShape)+
  theme_article()+geom_hline(yintercept = 0)
fe3.plot <- fe3.plot + geom_point(size=4, position=pos3) + 
  geom_errorbar(limits, width=0, colour="black", position = pos3) + 
  labs(x=element_blank(),y="Proportion") + 
  theme(legend.title=element_blank(),axis.text=element_text(size=10),
        legend.position = "bottom")

fe3.plot # Figure E.3
# Export Figure E.3
ggsave(file = "fe3.eps", width = 6, height = 3, dpi = 666)


